Simulation d'une cavité carrée avec les plans de symétrie (XZ) et (YZ), la postion du carré : Rotation 45 degré

In [1]:
import meep as mp
import matplotlib.pyplot as plt
import numpy as np
import copy

mp.quiet(True)
In [2]:
# Initialized parameter 

b = 1          # base of the square (On peut prend l'unité) 10 pixel 
pml = 0.02     # the thickness of pml layer
pad = 0.03     # the padding between the cavity and the pml layer
n = 1.5        # refractive index the cavity (used resin)
resolution = 200
len_scale = 25

sxy = b*np.sqrt(2) + 2*pml + 2*pad
cell_size = mp.Vector3(sxy, sxy, 0)

# Booudary conditions (using perfect matched layer)
boudary_layer = mp.PML(pml)

# Geometry of the cavity (rotated 45 degree)
geometry = [mp.Block(mp.Vector3(b,b,mp.inf),
                     mp.Vector3(1,1,0), 
                     mp.Vector3(1,-1,0),
                     center=mp.Vector3(),
                     material=mp.Medium(index=n))]

# Injected sources
wl = 4                       # wavelenght in micrometer / in vaccum        
fcen = len_scale/wl          # pulse center frequency
df = 20                      # pulse frequency width

pos_x = 0.25
pos_y = 0.41
src_centers = [mp.Vector3(-pos_x,pos_y), mp.Vector3(pos_x,pos_y), 
               mp.Vector3(pos_x,-pos_y), mp.Vector3(-pos_x,-pos_y)]
src_size = mp.Vector3(0)
src_amps = [1,1,1,1]
src_compt = mp.Ez

src = [mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[0], src_size, src_amps[0]),
       mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[1], src_size, src_amps[1]),
       mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[2], src_size, src_amps[2]),
       mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[3], src_size, src_amps[3])]
#Simulation
sim = mp.Simulation(cell_size= cell_size,
                    geometry= geometry,
                    sources=src,
                    symmetries=[mp.Mirror(mp.X), mp.Mirror(mp.Y)],
                    resolution=resolution,                    
                    boundary_layers=[boudary_layer])
In [3]:
# visualize the cavity
f = plt.figure(dpi=150)
sim.plot2D(ax = f.gca())
x0 = [-0.37]
y0 = [0.12]
plt.plot(x0, y0, "o")
plt.show()
In [4]:
# define step function to collect field Ez
def append_fields(sim):
    box = mp.Volume(center=pt_observed, size=mp.Vector3(0,0))
    ez_fields.append(sim.get_array(vol=box, component=mp.Ez))
In [5]:
ez_fields = []
pt_observed = mp.Vector3(-0.37,0.12)

sim.reset_meep()

sim.run(mp.after_sources(append_fields),
        until_after_sources=250)
In [7]:
ez0 = ez_fields.copy()
plt.figure()
plt.plot(ez0)
plt.title('Ez')
plt.show()
In [9]:
np.savez('data/Ez_even_t250', np.array(ez0))

Profil de modes de résonance

In [51]:
%matplotlib notebook
path = "/home/hong/Desktop/Microlaser1/data/"
file = "modes_Ez_even.npz"
npzfile = np.load(path + file)
modes = npzfile['arr_0']

plt.figure()
plt.scatter(modes[:,1], -modes[:,1]/2/modes[:,0])

plt.xlabel('$2\pi \omega_n$')
plt.ylabel('$Im(-2 \pi \omega_n)$')
plt.grid()
plt.title('modes de résonance')
x1 = [1.40173154e+01]
y1 = [-1.40173154e+01/2/3.08283178e+03]

x2 = [1.40185572e+01]
y2 = [-1.40185572e+01/2/6.11559092e+02]

x3 = [1.40256102e+01]
y3 = [-1.40256102e+01/2/2.09006114e+02]
plt.plot(x1, y1, "rs", x2, y2, "gs", x3,y3,"ys")
plt.show()
In [29]:
print(modes[(modes[:,1] > 14.5) 
            & (modes[:,1] < 15) 
            #& (modes[:,0]>100)
           ])

[[ 6.00977161e+01  1.49251142e+01  9.29913192e-06 -2.01761672e+00]
 [ 6.36095457e+01  1.47721797e+01  1.37801485e-05 -3.09161764e+00]
 [ 6.62814661e+01  1.46316631e+01  1.04391130e-05  2.69395338e+00]
 [ 7.01251239e+01  1.45000525e+01  1.94101271e-05 -1.26825982e+00]
 [ 7.74906983e+01  1.49410511e+01  6.26511345e-06  8.63106297e-01]
 [ 8.19920378e+01  1.48217736e+01  3.54381300e-05  1.57484938e+00]
 [ 8.94195536e+01  1.47212802e+01  6.35484726e-05  2.59834215e+00]
 [ 1.00672903e+02  1.46356352e+01  1.51437508e-04  1.03189754e+00]
 [ 1.17158205e+02  1.45642685e+01  1.15693130e-05 -1.32374701e+00]
 [ 1.47494085e+02  1.45050789e+01  3.73433625e-04  2.19154322e+00]
 [ 1.59100436e+02  1.49479460e+01  8.56727962e-05  2.07610684e+00]
 [ 2.67898518e+02  1.49186188e+01  3.88255690e-04  2.41258189e+00]
 [ 7.64024040e+02  1.49142306e+01  1.41252914e-03  2.88592162e+00]
 [ 3.14822797e+03  1.49120161e+01  8.91281767e-04 -1.98752382e+00]]

Méthode pour répéter

In [25]:
def run_mode(fcen, df, until) :
    # overwrite parameter of simulation
    fcen = fcen
    df = df
    src = [mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[0], src_size, src_amps[0]),
           mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[1], src_size, src_amps[1]),
           mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[2], src_size, src_amps[2]),
           mp.Source(mp.GaussianSource(fcen, fwidth=df), src_compt, src_centers[3], src_size, src_amps[3])]

    sim.reset_meep()       # reset everything
    sim.sources = src

    sim.run(mp.after_time(0, append_fields),
            until_after_sources = until)  # run until steady state

1. freq = 1.40173154e+01, famille #1

In [26]:
%matplotlib inline
In [37]:
pt_observed = mp.Vector3(-0.28,0.33)  # for function append_fields
In [28]:
ez_fields = []  # reset field
run_mode(fcen = 1.40173154e+01, df = 0.001, until = 100)
In [30]:
# Plot propagating field
ez1 = ez_fields.copy()
plt.figure()
#plt.axis([0.5,0.6,-17,17])
plt.plot(ez1)
plt.title('Ez')
plt.show()
In [35]:
# Plot field Ez in the last instant
plt.figure(figsize=(12,9))
x0 = [-0.28]
y0 = [0.33]
plt.plot(x0, y0, "s")
sim.plot2D(fields = mp.Ez)
Out[35]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fe2b6e92610>
In [36]:
freq = 1.40173154e+01
f = plt.figure(dpi=150)
Animate = mp.Animate2D(sim, fields=mp.Ez, f=f, realtime=False, normalize=True)
sim.run(mp.at_every(1/freq/20,Animate),until=1/freq)
plt.close()
Animate.to_jshtml(5)
Out[36]:

2. freq = 1.40185572e+01, famille #2

In [43]:
ez_fields = []
run_mode(fcen = 1.40185572e+01, df = 0.001, until = 100)
In [44]:
ez2 = ez_fields.copy()
plt.figure()
#plt.axis([0.5,0.6,-17,17])
plt.plot(ez2)
plt.title('Ez')
plt.show()
In [45]:
plt.figure(figsize=(12, 9))
x0 = [-0.28]
y0 = [0.33]
plt.plot(x0, y0, "s")
sim.plot2D(fields = mp.Ez)
Out[45]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fe2e851bc10>
In [46]:
freq = 1.40185572e+01
f = plt.figure(dpi=150)
Animate = mp.Animate2D(sim, fields=mp.Ez, f=f, realtime=False, normalize=True)
sim.run(mp.at_every(1/freq/20,Animate),until=1/freq)
plt.close()
Animate.to_jshtml(5)
Out[46]:

3. freq = 1.31178429e+01, famille #1

In [47]:
ez_fields = []
run_mode(fcen = 1.31178429e+01, df = 0.005, until = 200)
In [48]:
ez3 = ez_fields.copy()
plt.figure()
#plt.axis([0.5,0.6,-17,17])
plt.plot(ez3)
plt.title('Ez')
plt.show()
In [49]:
plt.figure(figsize=(12, 9))
x0 = [-0.28]
y0 = [0.32]
plt.plot(x0, y0, "s")
sim.plot2D(fields = mp.Ez)
Out[49]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fe268c5d6d0>
In [50]:
freq = 1.31178429e+01
f = plt.figure(dpi=150)
Animate = mp.Animate2D(sim, fields=mp.Ez, f=f, realtime=False, normalize=True)
sim.run(mp.at_every(1/freq/20,Animate),until=1/freq)
plt.close()
Animate.to_jshtml(5)
Out[50]: